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ABSTRACT 

A three-dimensional code for rotating blade-row flow analysis has been 
developed. The space discretization uses a cell-centered scheme with eigenvalues 
scaling for the artificial dissipation. The computational efficiency of a four-stage 
Runge-Kutta scheme is enhanced by using variable coefficients, implicit residual 
smoothing, and a full-multigrid method. 

An application is presented for the NASA rotor 67 transonic fan. Due to the 
blade stagger and twist, a zonal, non-periodic H-type grid is used to minimize the 
mesh skewness. The calculation is validated by comparing it with experiments in the 
range from the maximum flow rate to a near-stall condition. A detailed study of the 
flow structure near peak efficiency and near stall is presented by means of pressure 
distribution and particle traces inside boundary layers. 


* Assistant Professor 




INTRODUCTION 

In the last decade Fluid Dynamics has undergone impressive evolution both in the 
understanding and in the simulation of flow features. In this process. Computational 
Fluid Dynamics (CFD) is playing a more and more important role. Modem 
turbomachinery operates under very complex three-dimensional flow conditions, and 
further improvement requires detailed knowledge of the flow structure. Particularly, 
the need to estimate off-design conditions, secondary flows, and heat transfer forces 
us to look at viscous models. The real flow inside a turbomachine is unsteady and 
strongly influenced by rotor-stator interactions, wake, and clearance effects. Even if 
some important steps have been made in the time-accurate and time-averaged 
simulation of entire stages (e g. Rai, 1987, Rao and Delaney, 1990, Adamczyk et al., 
1990), a steady blade-row analysis should still be considered a basic tool for modem 
design. The works of Subramanian and Bozzola (1987), Chima and Yokota (1988), 
Choi and Knight (1988), Davis et al. (1988), Hah (1989), Nakahashi et al. (1989), 
Weber and Delaney (1991), and Dawes (1991) are some important steps in the 
prediction of three-dimensional viscous cascade flows. 

In 1988, the University of Florence started a joint project with NASA (ICASE 
and ICOMP) on viscous cascade flow simulation. During this research project, the 
TRAF2D and TRAF3D codes (TRAnsonic Flow 2D/3D) were developed (Amone et 
al. 1988, 1991, 1992). Those codes are capable of solving viscous cascade flows 
using H-type or C-type grids and of predicting heat transfer effects. In the present 
work, the procedure was extended to the case of rotating blade passages. Particular 
attention has been dedicated to important aspects such as minimization of the grid 
skewness, accuracy, and computational cost. 

As for accuracy, non-periodic C- or H-type grids are stacked in three dimensions. 
The removal of periodicity allows the grid to be only slightly distorted even for 
cascades having a large camber or a high stagger angle and twist. This allows us to 
pick up details of the throat flow with a reasonable number of grid points. A very low 
level of artificial dissipation is guaranteed by eigenvalues scaling, which is a three- 
dimensional extension of the one proposed by Swanson and Turkel (1987), and 
Martinelli and Jameson (1988). 

The two-layer eddy-viscosity model of Baldwin and Lomax (1978) is used for the 
turbulence closure. 

As for efficiency, the Reynolds-averaged Navier-Stokes equations are solved 
using a Runge-Kutta scheme in conjunction with accelerating techniques. Variable- 
coefficient implicit residual smoothing, as well as the Full- Approximation- Storage 
multigrid scheme of Brandt (1979) and Jameson (1983) are used in the TRAF3D 
code. Those accelerating strategies are implemented in conjunction with grid 
refinement to get a Full Multigrid Method. 

The code is validated by studying the NASA rotor 67 transonic fan. This fan has 
important viscous and three-dimensional effects and experiments are available in a 
wide range of flow rates from the maximum one to near stall. Moreover, this 
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geometry was recently proposed as an AGARD test case for code validation (Fottner, 
1990) and calculations from several authors are available for discussion. 

By using the accelerating strategies, accurate, viscous 3D solutions can be 
obtained in about half an hour on a modem supercomputer such as a Cray Y-MP. 


GOVERNING EQUATIONS 

Let p, u, v, w, p, T, E, and H denote respectively density, the absolute velocity 

components in the x, y, and z Cartesian directions, pressure, temperature, specific 

total energy, and specific total enthalpy. The three-dimensional, unsteady, Reynolds- 
averaged Navier- Stokes equations can be written for a rotating blade passage in 
conservative form in a curvilinear coordinate system £ jj, £ as, 

KJ' 1 Q) £F dG dH JF, , <?Gv ,dH* ;J (1) 

at at ap a$ a$ a tj a$ 

where the Cartesian system xy,z is rotating with angular velocity Q around the x 
axis. 
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The contravariant velocity components of eqs. (2) are written as. 
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V = T}, + T] x u + TJ y V +J] 2 w 

We'Z+C'U+CyVK W 


( 3 ) 


and the transformation metrics are defined by, 

4L =J(y n zs-y ( z n ), 4 y =J (z n x{-z4x n ) 

4 z =J(x n y ( - x 4 y„) , Vx =J ( y 4 z 4 -y&) 

%=J( z 4 x 4 ~ z 4 x 4 ), r] x =J( Xc y 4 -x 4 y 4 ) 

C =J(y 4 Zr,-Z 4 y„), Gy=J(z4X n -X 4 zJ (4) 

Cz ~ X 4 y JJ ~ y 4 Xrj) J y,Zy Zt^z 

n, = -x, Vx -y, % -z, rjz, £ = -X'C x -yXy -Z'Cz 
x <=0, y t =-Oz, Zt=&y 

where the Jacobian of the transformation J is, 

T 1 =x 4 y,,z 4 +x n y ( z 4 +x 4 y 4 zn-x 4 y 4 z v -x n y 4 z 4 -x 4 y n z 4 (5) 

The viscous flux terms are assembled in the form, 
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Txx = 2 MUx+ 2,{ux +Vy+ w> z ) 
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T a = 2H\Vi +2.(lix + Vy 
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Px - U Txx + V Txy + W Txz +k Tx 
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P 2 =U Tzx +V Tzy +W Tzz +kTz 


and the Cartesian derivatives of (7) are expressed in terms of tj-, and ^-derivatives 
using the chain rule, i.e.. 


Ux=£, x Ut + r\xUn+$ x uz 


( 8 ) 
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The pressure is obtained from the equation of state, 

p =pRT (9) 

According to the Stokes hypothesis, X is taken to be -2 p/3 and a power law is used to 
determine the molecular coefficient of viscosity p as function of temperature. The 
eddy- viscosity hypothesis is used to account for the effect of turbulence. The 
molecular viscosity p and the molecular thermal conductivity k are replaced with. 


P=Pi+P, 0 °) 



where c p is the specific heat at constant pressure, Pr is the Prandtl number, and the 
subscripts / and / refer to laminar and turbulent respectively. The turbulent quantities 
p t and Pr t are computed using the two-layer mixing length model of Baldwin and 
Lomax (1978). The contribution of the eddy viscosity is computed separately in the 
blade-to-blade direction 7 and in the spanwise direction £ The inverse of the square 
of the wall distances d is then used to compute the resulting eddy viscosity, 


1 


f = 



( 12 ) 


0 3 ) 

The transitional criteria of Baldwin and Lomax is adopted on the airfoil surface while 
on the end walls, the shear layer is assumed to be fully-turbulent from the inlet 
boundary. 


SPATIAL DISCRETIZATION 

Traditionally, using a finite-volume approach, the governing equations are 
discretized in space starting from an integral formulation and without any intermediate 
mapping (e.g. Jameson at al., 1981, Ni, 1981, Holmes and Tong, 1985). The 
transformation metrics of (4) can be then associated with the projections of the face 
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areas as the contravariant components of (3) can be related to the normal components 
of the relative velocity. In the present work, due to the large use of eigenvalues and 
curvilinear quantities, we found it more convenient to map the Cartesian space ( x,y,z ) 
in a generalized curvilinear one (4, rj, 0. In the curvilinear system, the equation of 
motion (1) can be easily rewritten in integral form by means of Green's theorem and 
the metric terms are handled following the standard finite-volume formulation. The 
computational domain is divided into hexahedrons and the transformation metrics are 
evaluated so that the projected areas of the cell-faces are given by the ratio of the 
appropriate metric derivatives to the Jacobian ones, i.e. &J is the projection onto the 
x-axis of a cell face at a fixed 4 location. A cell-centered scheme is used to store the 
flow variables. On each cell face the convective and diffusive fluxes are calculated 
after computing the necessary flow quantities at the face center. Those quantities are 
obtained by a simple averaging of adjacent cell-center values of the dependent 
variables. 


BOUNDARY CONDITIONS 

In cascade calculations we have four different types of boundaries: inlet, outlet, 
solid walls, and periodicity. At the inlet, the presence of boundary layers, on hub and 
tip end walls, is accounted for by giving a total pressure and a total temperature 
profile whose distribution simulates the experimental one. According to the theory of 
characteristics, the flow angles, total pressure, total temperature, and isentropic 
relations are used at the subsonic-axial inlet, while the outgoing Riemann invariant is 
taken from the interior. At the subsonic-axial outlet, the average value of the static 
pressure at the hub is prescribed and the density and components of velocity are 
extrapolated together with the circumferential distribution of pressure. The radial 
equilibrium equation is used to determine the spanwise distribution of the static 
pressure. On the solid walls, the pressure is extrapolated from the interior points, and 
the no-slip condition and the temperature condition are used to compute density and 
total energy. For the calculations presented in this paper, all the walls have been 
assumed to be at a constant temperature equal to the total inlet one. 

Cell-centered schemes are generally implemented using phantom cells to handle 
the boundaries. The periodicity from blade passage to blade passage is, therefore, 
easily overimposed by setting periodic phantom cell values. On the boundaries where 
the grid is not periodic, the phantom cells overlap the real ones. Linear interpolations 
are then used to compute the value of the dependent variables in phantom cells. 
Even if this approach does not guarantee conservation of mass, momentum, and 
energy, no accuracy losses have been experienced unless strong flow gradients occur 
along non-periodic grid boundaries with strong differences in mesh size. In the 
present work, the number of fine cells to be used for interpolation on a coarse cell 
never exceeded three. 
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The clearance region is handled by imposing periodicity conditions across the 
airfoil without any modellization of the blade cross-section. 


ARTIFICIAL DISSIPATION 

In viscous calculations, dissipating properties are present due to diffusive terms. 
Away from the shear layer regions, the physical diffusion is generally not sufficient to 
prevent the odd-even point decoupling of centered schemes. Thus, to maintain 
stability and to prevent oscillations near shocks or stagnation points, artificial 
dissipation terms are also included in the viscous calculations. Equation (1) is written 
in semi-discrete form as. 


-%+C(Q)-D(Q)=0 


(14) 


where the discrete operator C accounts for the physical convective and diffusive 
terms, while D is the operator for the artificial dissipation. The artificial dissipation 
model used in this paper is basically the one originally introduced by Jameson, 
Schmidt, and Turkel (1981). In order to minimize the amount of artificial diffusion 
inside the shear layer, the eigenvalues scaling of Martinelli and Jameson (1988), and 
Swanson and Turkel (1987) have been used to weight these terms. The quantity 
D(Q) in eq. (14) is defined as. 


D(Q) =($ ~D\ +D 2 r,-D 4 r, +$ ~D$Q 


where, for example, in the £ curvilinear coordinates we have, 

D\0 = V ? U i+]/2,j,k dH/2.].k) Al,Q ljk 
D 4 sQ=V'U +l/2J.k d+i/2.j.k) A4 VfA^Q I jk 


(15) 


(16) 


i,j,k are indices associated with the £, 77 , 4 " directions and V$, At, are forward and 
backward difference operators in the £ direction. The variable scaling factor A is 
defined for the three-dimensional case as, 

Am/2j* ( 17 ) 

where, 

A4 = 0^ (l 8 ) 
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The definition of the coefficient <2>has been extended to the three-dimensional case as 
follows, 


0 4 =1 + 
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where X^, and X^ are the scaled spectral radii of the flux Jacobian matrices for the 
convective terms. 


2,=| f I -mjti’ + rt’ + rf (20) 

and a is the speed of sound. Note that the effect of the grid rotation is accounted for 
in (20) through the definition of the contravariant components of velocities of (3). 
The exponent a is generally defined by 0<o <1 , and for two-dimensional 
applications, a value of 2/3 gives satisfactory results. In three-dimensional cascade 
flow calculations, we generally have highly stretched meshes in two directions near 
comers. We found that <7=0.4 introduces enough scaling without compromising the 
robustness. The coefficients and d 4) use the pressure as a sensor for shocks and 
stagnation points, and are defined as follows. 


£j /Zhk =k (2 >MAX( Vi-},j,k > VS ,j,k* Vi+l,j,k> Vi+2,j,k) 


( 21 ) 
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( 22 ) 
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43aj* - MAX\o,{k w - 


(23) 


where typical values for the constants KP> and K< 4 > are 1/2 and 1/64 respectively. For 
the remaining directions 77 and £ the contribution of dissipation is defined in a similar 
way. The computation of the dissipating terms is carried out in each coordinate 
direction as the difference between first and third difference operators. Those 
operators are set to zero on solid walls in order to reduce the global error on the 
conservation property and to prevent the presence of undamped modes (Pulliam, 
1986, and Swanson and Turkel, 1988,). 

It is important to anticipate now that from the definition of residual of (25), 
variable scaling, and time steps of (26,27,28), it results that the artificial dissipation is 
scaled with a factor proportional to the ratio between the global time step and the 
inviscid time step. Close to solid walls, the grid volume is very small and viscous time 
step limitation is dominant. The ratio of the time step over the inviscid one becomes 
very small and most of the artificial dissipation is removed. 


TIME-STEPPING SCHEME 

The system of the differential equation (14) is advanced in time using an explicit 
four stage Runge-Kutta scheme until the steady-state solution is reached . A hybrid 
scheme is implemented, where, for economy, the viscous terms are evaluated only at 
the first stage and then frozen for the remaining stages. If n is the index associated 
with time we will write it in the form. 


Qf 0> =Q n 

Qf" =<2 0) +aiA$ 0) ) 

a 7 *-#* +<*!{&>) 

$ 3) =<2 0> +a 3 A<2 2) ) ( 24 ) 

<2 4) =<2 0) +cuA&) 

Q n+ ' =& 4> 


where the residual R (Q) is defined by, 

R(Q)=AtJ[C(Q)-D(Oj\ (25) 
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Good, high-frequency damping properties, important for the multigrid process, have 
been obtained by performing two evaluations of the artificial dissipating terms, at the 
first and second stages. It is worthwhile to notice that, in the Runge-Kutta time- 
stepping schemes, the steady state solution is independent of the time step; therefore, 
this stepping is particularly amenable to convergence acceleration techniques. 


ACCELERATION TECHNIQUES 

In order to reduce the computational cost, four techniques are employed to speed 
up convergence to the steady state-solution. These techniques: 1) local time- 
stepping; 2) residual smoothing; 3) multigrid; 4) grid refinement; are separately 
described in the following. 


Local Time-Stepping 

For steady state calculations with a time-marching approach, a faster expulsion of 
disturbances can be achieved by locally using the maximum available time step. In the 
present work the local time step limit At is computed accounting for both the 
convective (AtJ and diffusive (At 4 ) contributions as follows. 


At =co 

V 


At+At d ) 


(26) 


where c 0 is a constant usually taken to be the Courant-Friedrichs-Lewy (CFL) 
number. Specifically, for the inviscid and viscous time step we used, 

At c = (27) 


At a = (28) 

K.-^ASnSf +S ! ( S( +S(S$ 

prr 

where / is the specific heat ratio and, 

S\=x\+y\+z\, S 1 r t = x 2 n +y n +z\, s\=x\ + y\+z\, (29) 
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Kf being a constant whose value has been set equal to 2.5 based on numerical 
experiments. 


Residual Smoothing 

An implicit smoothing of residuals is used to extend the stability limit and the 
robustness of the basic scheme. This technique was first introduced by Lerat in 1979 
in conjunction with Lax-Wendroff type schemes. Later, in 1983, Jameson 
implemented it on the Runge-Kutta stepping scheme. In three dimensions we carried 
out the residual smoothing in the form, 

(l ~P ( VtA&l - P , V„A&1 ~Pi V^R =R (30) 

where the residual R includes the contribution of the variable time step and is defined 

by (25) and R~ is the residual after a sequence of smoothing in the £ 7 , and £ 
directions with coefficients 0$ , 0^ , and 0$. For viscous calculations on highly 
stretched meshes the variable coefficient formulations of Martinelli and Jameson 
(1988) and Swanson and Turkel (1987) have proven to be robust and reliable. In the 
present paper, the expression for the variable coefficients 0 of (30) has been modified 
to be used in three dimensions as follows, 


0 c =MAX{ 
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where the coefficients 0 r 0 y and 0 Z are the ones defined in eqs. (19), and CFL, and 
CFL * are the Courant numbers of the smoothed and unsmoothed scheme 

respectively. For the hybrid four-stage scheme we used CFL=5, and CFL‘=2.5. 
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Multiqrid 

This technique was developed in the beginning of the 1970s for the solution of 
elliptic problems (Brandt, 1979) and later was extended to time-dependent 
formulations (Ni, 1981, and Jameson, 1983). The basic idea is to introduce a 
sequence of coarser grids and to use them to speed up the propagation of the fine grid 
corrections, resulting in a faster expulsion of disturbances. In this work we used the 
Full Approximation Storage (FAS) schemes of Brandt (1979) and Jameson (1983). 

Coarser, auxiliary meshes are obtained by doubling the mesh spacing and the 
solution is defined on them using a rule which conserves mass, momentum, and 
energy, 


( 32 ) 

where the subscripts refer to the grid spacing, and the sum is over the eight cells 
which compose the 2h grid cell. Note that this definition coincides with the one used 
by Jameson when the reciprocal of the Jacobians are replaced with the cell volumes. 
To respect the fine grid approximation, forcing functions P are defined on the coarser 
grids and added to the governing equations. So, after the initialization of Q, 2 h using 
eq.(32), forcing functions Pjh are defined as, 


P2h = ERtfQh) ~ P™(@2h) ( 33 ) 

and added to the residuals fch to obtain the value which is then used for the 
stepping scheme. 


Fbk=R2H(Q2 h )+P2H ( 34 ) 

This procedure is repeated on a succession of coarser grids and the corrections 
computed on each coarse grid are transferred back to the finer one by bilinear 
interpolations. 

A V-type cycle with subiterations is used as a multigrid strategy. The process is 
advanced from the fine grid to the coarser one without any intermediate interpolation, 
and when the coarser grid is reached, corrections are passed back. One Runge-Kutta 
step is performed on the h grid, two on the 2h grid, and three on all the coarser grids. 
It is our experience in cascade flow calculations that subiterations increase the 
robustness of the multigrid. 

For viscous flows with very low Reynolds number or strong separation, it is 
important to compute the viscous terms on the coarse grids, too. The turbulent 
viscosity is evaluated only on the finest gnd level and then interpolated on coarse 
grids. 
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On each grid, the boundary conditions are treated in the same way and updated at 
every Runge-Kutta stage. For economy, the artificial dissipation model is replaced on 
the coarse grids with constant coefficient second-order differences. 

The interpolations of the corrections introduce high frequency errors. In order to 
prevent those errors from being reflected in the eddy viscosity, turbulent quantities are 
updated after performing the stepping on the fine grid. On coarse grids, the turbulent 
viscosity is evaluated by averaging the surrounding fine grid values. 

Grid Refinement 

A grid refinement strategy is used to provide a cost-effective initialization of the 
fine grid solution. This strategy is implemented in conjunction with multigrid to 
obtain a Full Multigrid (FMG) procedure. With the FMG method, the solution is 
initialized on a coarser grid of the basic grid sequence and iterated a prescribed 
number of cycles of the FAS scheme. The solution is then passed, by bilinear 
interpolations, onto the next, finer grid and the process is repeated until the finest grid 
level is reached. In the present paper we have introduced three levels of refinement 
with respectively two, three, and four grids. 


COMPUTATIONAL GRID 

The three-dimensional grids are obtained by stacking two dimensional grids 
generated on blade-to-blade surfaces at constant radii (£7 plane). In order to 
minimize the grid skewness, in the blade-to-blade projection, the grids structure can 
be chosen on the basis of the blade geometry and flow conditions. Turbine blades are 
generally characterized by blunt leading edge and high turn with a subsonic incoming 
flow in the relative plane. For these geometries, non-periodic C-type grids, have 
proven to be effective (e.g. Amone et al., 1991, 1992). In the case of a compressor, 
and particularly for a fan, the leading aspect of the geometry is the high stagger and 
twist, while the leading edge is often quite sharp (fig. 1). In addition, the incoming 
flow in the relative plane can be supersonic for a large part of the blade span. As 
consequence, a C-type structure of the grid may smear too much of the bow shock 
away from the leading edge (i.e. Weber and Delaney, 1991). In the present work, a 
non-periodic H-type grid was implemented. The removal of mesh periodicity allow 
the grid to accommodate highly staggered airfoils with a low level of skewness. To 
minimize the undesired interaction between strong flow gradients and non-periodic 
boundaries the mesh correspondence is broken before the leading edge and on the 
wake, but not inside the blade channel. The inviscid grids are elliptically generated, 
controlling the grid spacing and orientation at the wall. Viscous blade-to-blade grids 
are then obtained from inviscid grids by adding lines near the wall with the desired 
spacing distribution. 
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For highly twisted blades, a low grid skewness in the various spanwise sections 
can also be maintained by adjusting the grid-point distributions on the suction and 
pressure sides of the blade. 

In the spanwise direction (£) a standard H-type structure is used. Near the hub 
and tip endwalls, geometric stretching is used for a specified number of grid points, 
after which the spanwise spacing remains constant. 


APPLICATION AND DISCUSSIONS 

As validation of the procedure that has been described above, the TRAF3D code 
was used to study the NASA rotor 67 transonic fan. This first-stage rotor of a two- 
stage transonic fan was designed and tested with laser anemometer measurements at 
NASA Lewis (Pierzga and Wood, 1985). The rotor has 22 low aspect ratio ( 1.56) 
blades and was designed for a rotational speed of 16043 rpm, with a total pressure 
ratio of 1.63 and a mass flow of 33.25 kg/s. Experiments for simple blade-row code 
validation are available for the rotor without inlet guide vanes or downstream stators, 
and include detailed data near the peak efficiency and stall conditions. This rotor 
geometry has been recently proposed as an AGARD test case (Fottner, 1990) and 
several authors (Adamczyk et al. 1991, Chima, 1991, Hah and Reid, 1992, Jennions 
and Turner, 1992) have computed this geometry trying to understand the complex 
nature of transonic rotor flows. Therefore theoretical predictions from different codes 
are also available in the bibliography for comparisons and discussions. 

A three-dimensional view of an inviscid grid for the rotor is given in fig. 1 . In 
previous works by the author (Amone et al., 1991, 1992), a two- and three- 
dimensional grid dependency study was earned out in order to figure out the mesh 
requirements necessary to obtain a space-converged calculation, especially for viscous 
details such as losses, skin friction and heat transfer. Those results can be 
extrapolated to the case of rotating blade passages. Using the algebraic turbulence 
model of Baldwin and Lomax (1978), a y + at the wall of about four, with a Reynolds 
number of about one million, gives satisfactory results unless heat transfer details are 
needed. This y + value was achieved with a mesh spacing at the wall in the blade-to- 
blade direction of 2x1 O' 4 times the hub axial-chord. In the spanwise direction, due to 
the relatively thick inlet boundary layer, the mesh spacing at the wall was fixed at 
1*10~ 3 times the hub axial-chord. One hundred thirty-seven points were used in the 
streamwise direction and 49 in the blade-to-blade and hub-to-tip directions. Sixty-five 
points were located on the suction and pressure side of the airfoil. In the spanwise 
direction, four cells lie inside the clearance region. Three grid sections at 70%, 30%, 
and 10%> of the span from the shroud, are shown in fig. 2(b), (c), and (d), 
respectively. Figure 2 (a) gives a meridional view of the grid. 

Following the approach suggested by Pierzga and Wood (1985) the comparison 
between calculations and experiments is carried out using the mass flow rate 
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nondimensionalized with the choke flow rate as equivalence criteria. However, good 
agreement was also found in terms of absolute quantities. The TRAF3D predicts 34.5 
kg/s while 34.96 kg/s was the measured one. This underestimation of the choke flow 
rate of about 1.3% agrees with the viscous prediction of Chima (1991) and Jennions 
and Turner (1992). The peak efficiency and near stall conditions correspond to a 
nondimensional flow rate of . 989 and .924, respectively. 

The inlet boundary layer in the endwall region is accounted for by giving a total 
pressure profile. The thickness of the boundary layer is taken from experiments and 
the 1/7 power law velocity profile is used to estimate the distribution of total pressure. 
The core flow is assumed uniform. 

The convergence of the root mean square of the norm of residuals is given in fig. 
3. This calculation refers to the near-peak efficiency condition, and requires about 35 
minutes on the NASA Lewis Cray Y-MP to achieve a four decades' reduction in the 
residuals and corresponds to 100 multigrid cycles on the finest grid level using four 
grids. More than 35 minutes were needed only for solutions close to the stall 
condition. The global mass flow error has always been found to be less then 10r$ . 

Calculations were performed for 13 different values of the nondimensional flow 
rate in order reproduce the operating characteristics of the rotor at the design speed. 
Results are summarized in fig. 4 in terms of rotor adiabatic efficiency and rotor total 
pressure ratio. A nondimensional value of the mass flow rate equal to about .91 was 
the smaller value for which a steady solution was obtained. Further increase in the exit 
pressure would produce tip stall. The prediction of the rotor efficiency is quite good 
and agrees with the indication of Jennions and Turner (1992) in the fact that the 
abrupt decrease going from peak efficiency to stall conditions is smoothed out in the 
calculations. The peak efficiency is predicted at a slightly lower mass flow rate. 
Some underestimation in the total pressure ratio across the rotor should be noticed, 
but it seems to be evident only when approaching the stall condition. It is believed 
that the following aspects could be investigated in order understand and/or address 
this problem: 

• the distribution of the inlet boundary layer has been modelled only in terms of 
boundary layer thickness while experiments also show some gradients in the 
inviscid core of the inlet flow. 

• the Baldwin-Lomax turbulence model is not very effective for transonic flow 
with strong adverse pressure gradients (Dawes, 1990). 

• measurements of the machine geometry while running have shown some 
deformation not included in the present calculations (Fottner, 1990). 

• the tip clearance model which has been used is quite simple and clearance flow 
becomes important when the stall condition is approached (Adamczyk et al., 
1991, and Jennions and Turner, 1992). 


16 



The quality of the computed solutions is evaluated by comparing the spanwise 
distribution of circumferential energy-averaged thermodynamic quantities downstream 
of the rotor to experiments. The predicted distribution of static pressure, flow angle, 
total pressure, and total temperature at the rotor exit is compared to experiments in 
figs. 5 and 6 for conditions near peak efficiency and stall respectively. The agreement 
is good on the whole. The radial distribution of the static pressure is well reproduced 
(see figs. 5 (a) and 6 (a)), which indicates a good estimation of the global losses. The 
exit angle is up to five degrees off in the central part of the blade span (figs. 5 (b) and 
6 (b)), but this agrees with the calculations of Chima (1991) and Jennions and Turner 
(1992). From the plots of total pressure and total temperature of figs. 5 and 6 (c) and 
(d) we can see that the overall prediction of the rotor characteristic is well 
reproduced. 

As well known, the flow pattern inside a rotor like the NASA 67 transonic fan is 
very complicated and characterized by effects such as shock-boundary layer 
interaction, clearance flow, and three-dimensional separation with vortices roll up. 
One attempt at interpreting the structure of the computed flow field can be analyzing 
the relative flow on blade-to-blade and meridional surfaces. Pictures of the limiting 
streamlines inside the boundary layer can be obtained by means of particle traces. 

Blade-to-blade flow 

Figures 7 and 8 report the computed and measured relative Mach number 
contours at three different locations of the blade span. The agreement with 
experiments is qualitatively good and the bow shock is not too smeared away from 
the leading edge. It is also interesting to notice that those contours agree very well 
with the ones obtained by Jennions and Turner (1991) using a two-equation model for 
the turbulence closure. 

Near peak efficiency, the shock system has a lambda structure with a bow shock. 
The passage shock crosses the blade channel and involves about 30% of the upper 
part of the airfoil (see figs. 7, and also 10, and 11). Approaching the stall condition 
(fig. 8) the passage shock moves upstream and stands in front of the blade so that the 
airfoil pressure side is no longer intercepted (see also figs. 12 and 13). 

Blade surface flow pattern 

Pressure distribution and a restriction of the particle traces close to the airfoil 
surface are used to interpret the flow pattern inside the blade boundary layer. A 
schematic of this structure for the peak efficiency condition is depicted in fig. 9. As 
discussed by Weber and Delaney (1991) and Hah and Reid (1991), most of the 
separation and outward flow is observed on the suction side of the blade. The 
passage shock is quite strong in the upper part of the rotor and the losses in axial 
momentum result in a rapid turn toward the shroud. The separation due to shock- 
boundary layer interaction is evident in figs. 10 and 11. Separation lines are 
characterized by flows going toward the line, while where the flow reattaches, the 
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particle traces look like they are going away from the line (see. fig. 9). Such a 
situation is very clear in figs. 10 (c) and 13 (d). The passage shock also induces 
separation on the blade pressure side but only in the very upper part of the airfoil. 
The blow up of fig. 10 (c) shows the abrupt radial migration with separation and 
reattachment. 

In the central part of the blade span, the passage shock loses intensity and the 
flow lift off is mostly related to the adverse pressure gradient in the axial direction. 
Eventually, near the trailing edge, the flow separates on the suction side (see fig. 
11(b)). 

Close to the blade root, the flow is strongly influenced by a vortex roll up on the 
leading edge. As also pointed out by Chima (1991), particles undergo a high relative 
incidence close to the hub, which causes the low-momentum fluid to separate and 
migrate radially outward ( see figs. 10 (b) and 1 1 (b)). 

In accordance with the calculations of Weber and Delaney (1991), Chima (1991), 
and Jennions and Turner (1991), a separation bubble is observed on the blade suction 
side near the hub (figs. 9, 1 1 (b), and 13 (d)). 

The flow structure of the near-stall operating condition is qualitatively similar to 
the peak efficiency one previously discussed. Now the passage shock has moved 
upstream and the pressure side is shock free (fig. 12 (a)) with basically no radial flow 
mixing (fig. 12 (b)). On the contrary, on the suction side, the shock involves most of 
the blade span (fig. 13 (a)) and induces a strong outward flow with very clear 
separation and reattachment lines (fig. 13 (d)). The vortex roll up on the pressure 
side of fig. 10 (d) has now moved upstream and stands in front of the blade leading 
edge (fig. 13 (c)). The separation bubble on the suction side near the hub looks 
slightly smaller. On the shroud, the bow shock interacts with the casing boundary 
layer and the flow separates as depicted in fig. 1 3 (b). 

Clearance flow pattern 

A picture of the clearance flow pattern is obtained by plotting the relative Mach 
number contours and a restriction of the particle traces on a blade-to-blade surface 
midway between the blade tip and the casing. Once again the flow structure looks 
very similar to the one computed by Jennions and Turner (1991). 

The TRAF3D code predicts two tip vortices which intersect before mid channel. 
At peak efficiency (fig. 14 (a) and (b)), the shock system still shows a lambda 
structure and interacts with the tip vortices. A first vortex is observed close to the 
leading edge and a second leakage vortex forms at an axial location just after the 
pressure side shock. The two vortices sum up before crossing the passage shock as 
shown in fig. 14. When approaching the stall condition there is interaction between 
the leading edge shock separation and the leading edge vortex (Adamczyk at. al 
( 1 99 1 ), and Jennions and Turner (1991)). The leading edge vortex is now associated 
with the casing separation of the bow shock, while the leakage vortex has moved 
upstream. The two vortices interact very soon while going towards the pressure side 
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of the next consecutive blade (fig. 5(b) and (c). Figure 15 (d) indicates a strong link 
between the shock system and the tip vortices for this flow condition. 

Hub endwall flow pattern 

Particle traces in the hub boundary layer are reported for completeness in figs. 16 
and 17. The high angle of attack experienced by the flow in this region and 
previously discussed is evident. The separation bubble close to the trailing edge 
causes a easily visible lack of particles, which, when injected close to the blade, roll up 
radially. 


CONCLUSIONS 

The central-difference, finite-volume scheme with eigenvalues scaling for artificial 
dissipation terms, variable-coefficient implicit smoothing, and full multigrid has been 
extended to predict three-dimensional rotating blade passages. The procedure has 
been validate by comparing it with experiment for the NASA rotor 67 in a wide range 
of mass flow rate. With these accelerating strategies, detailed three-dimensional 
viscous solutions can be obtained for a reasonable fine grid in about half an hour on a 
modem supercomputer. 
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Fig. 1 Three-dimensional inviscid grid for the NASA rotor 67 
transonic fan. 
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c) grid at 30% span d) grid at 10% span 

from the shroud from the shroud 

Fig. 2 137x49x49 viscous grid for the NASA rotor 67 
transonic fan. 



Fig. 3 Convergence history for the near peak efficiency 
condition (NASA rotor 67 transonic fan). 


24 


























1.30 


1.25 


o 1-20 


< 1.151 

■ o . o o 

^ 1.10 

- 

1.05' 

1.00 

1 1 1 1 


0.0 0.2 0.4 0.6 0.8 1.0 


^hub)/ (^tip ^hub) 


Fig. 6 Predicted and measured exit survey data near stall 
(NASA rotor 67). 
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Fig. 7 Measured (left) and predicted (right) relative Mach 

number contours near peak efficiency (NASA rotor 67). 
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clearance flow 



Fig. 9 Schematic of the suction side boundary-layer flow 
structure near peak efficiency (NASA rotor 67). 


suction side shock 



Fig. 10 Pressure contours and particle traces close to the blade 
pressure side near peak efficiency (NASA rotor 67). 
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pressure side shock 



a) pressure b) particle traces 


Fig. 1 1 Pressure contours and particle traces close to the blade 
suction side near peak efficiency (NASA rotor 67). 


suction side shock 



a) pressure 



b) particle traces 


Fig. 12 Pressure contours and particle traces close to the blade 
pressure side near stall (NASA rotor 67). 
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a) pressure 


particle traces 

Fig. 13 Pressure contours and particle traces close to the blade 
suction side near stall (NASA rotor 67). 










Fig. 16 Particle traces close to the hub end wall 
near peak efficiency (NASA rotor 67). 



Fig. 17 Particle traces close to the hub endwall 
near stall (NASA rotor 67). 
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